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Abstract 


The straightforward automatic-differentiation and the hand-differentiated incremental iterative methods are 
interwoven to produce a hybrid scheme that captures some of the strengths of each strategy. With this compromise, 
discrete aerodynamic sensitivity derivatives are calculated with the efficient incremental iterative solution algorithm 
of the original flow code. Moreover, the principal advantage of automatic differentiation is retained (i.e., all 
complicated source code for the derivative calculations is constructed quickly with accuracy). The basic equations 
for second-order sensitivity derivatives are presented; four methods are compared. Each scheme requires that 
large systems are solved first for the first-order derivatives and, in all but one method, for the first-order adjoint 
variables. Of these latter three schemes, two require no solutions of large systems thereafter. For the other two 
for which additional systems are solved, the equations and solution procedures are analogous to those for the first- 
order derivatives. From a practical viewpoint, implementation of the second-order methods is feasible only with 
software tools such as automatic differentiation, because of the extreme complexity and large number of terms. 
First- and second-order sensitivities are calculated accurately for two airfoil problems, including a turbulent- 
flow example; both geometric-shape and flow-condition design variables are considered. Several methods are 
tested; results are compared on the basis of accuracy, computational time, and computer memory. For first-order 
derivatives, the hybrid incremental iterative scheme obtained with automatic differentiation is competitive with the 
best hand-differentiated method; for six independent variables, it is at least two to four times faster than central 
finite differences and requires only 60 % more memory than the original code; the performance is expected to 
improve further in the future. 
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1.0 Introduction 


The use of advanced computational fluid dynamics (CFD) analysis codes in multidisciplinary design opti- 
mization (MDO) studies and applications via sensitivity analysis requires the efficient and accurate calculation 
of individual-discipline sensitivity derivatives (SD). The incremental iterative method (IIM) was proposed and 
demonstrated to provide such first-order (FO) SD from a two-dimensional (2-D) thin-layer Navier-Stokes code 
(TLNS) for both geometric (shape) and nongeometric (flow) design variables in Refs. [1] and [2]. The IIM allows 
accurate, consistent discrete SD to be obtained with computational efficiency (with respect to both computational 
time and memory requirements). Furthermore, the DM also allows the use of approximate matrix operators for 
further efficiency, parallelization, or robustness, etc. Results for FO SD from an DM for three-dimensional (3- 
D) Euler codes (Refs. [3]-[5]) have also been presented. In all of the above cited works, the discretized flow 
residuals were differentiated by hand (also called the quasi-analytical (QA) method) and assembled to obtain the 
FO SD by an IIM. 

In the present study, numerical results are given for the application of automatic differentiation (AD) (Refs. 
[6]-[8]) to obtain FO aerodynamic SD from an IIM for the same 2-D TLNS code and sample problems studied 
in Ref. [1]. The numerical results are compared on the basis of accuracy and computational time and memory. 
Previous FO SD from hand-differentiated QA IIM and central finite-difference (CD) methods (Ref. [1]) are 
compared with newly obtained AD results for both the straightforward and IIM form applications. This latter 
approach is new; previous straightforward applications of AD to advanced CFD codes (Refs. [9]-[12]) did not 
result in IIM forms, as will be discussed subsequently. This problem was recognized and noted in Refs. [2] and 
[9]-[l 1], in which the use of AD in IIM forms was proposed. 

An additional focus of this study is the development of the basic equations for computing second-order 
(SO) discrete aerodynamic SD, which yields four methods; where applicable, the incremental iterative forms of 
these equations are also given. Numerical results are shown for the same 2-D sample problems for which FO 
SD are calculated. 

The SO aerodynamic SD are of interest for several reasons. For example, aerodynamic stability derivatives are 
required by the controls discipline as an input; therefore, inclusion of controls in a gradient-based MDO procedure 
means that the sensitivities of stability derivatives are needed, which are SO SD. Secondly, in constructing function 
approximations for nonlinear flow behavior, the expansions that use FO derivative information are only of limited 
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usefulness. For N independent variables D, the truncated Taylor sales 
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di 


f(D+AD) ~ f(D) + Y, 

jtr dD > 
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is a linear approximation. If the derivatives 3 ^ 5 - are available, then 


df 


f(D+AD) ~ f(D) + ^— ADj 
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may exhibit much of the nonlinear behavior of f. Thirdly and more importantly, SO optimization techniques 
can be employed. 

The remainder of the paper is organized as follows: Section 2 discusses first-order derivatives, including 
development of the equations, methods, results (tabulated in Appendix A), and conclusions. Section 3 discusses 
second-order derivatives, including development of the equations, methods, results (tabulated in Appendix B), 
and conclusions. Section 4 presents a summary and final conclusions. For the convenience of the reader, the 
acronyms used throughout the paper are collected in Appendix C. 


2.0 First-Order Sensitivity Derivatives (FO SD) 

2.1 Basic Equations and Incremental Iterative Forms 

A brief review is given of the basic equations of FO discrete aerodynamic sensitivity analysis; also included 
is the incremental iterative form for solving the equations. A thorough discussion is given in Refs. [1] and [2]. 
Reference [13] provides an overview of recent advances in FO sensitivity analysis for CFD. 

After discretization, the nonlinear, multidimensional steady-state governing equations of fluid flow and the 
boundary conditions are approximated as a large system of coupled nonlinear algebraic equations as 

R = R(Q(b),X(b),b) = 0 (2.1) 

where Q is the vector of field variables, X is the computational grid, and b is the vector of independent input 
(design) variables. Similarly, the vector of aerodynamic output functions F is dependent on Q, X, and b as 

F = F(Q(b),X(b),b) (2.2) 

In Eqs. (2.1) and (2.2) and all subsequent equations, all applicable terms are evaluated at the steady-state flow 
conditions, unless explicitly superscripted with an appropriate iteration index. 
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2.1.1 The Direct Differentiation (DD) Method. Differentiation of Eqs. (2.2) and (2.1) with respect to b 
yields the respective matrix equations 


d F' = — Q ' + —X' + — 
9Q 9X 9b 

r'= 5Sq. + Sx' + |£=o 
ax ab 


(2-3) 

(2.4) 


where °F' = g; R'= ft; Q'= f ; and X'= &. 

The matrix D F' contains the sensitivity derivatives of interest; the superscript D denotes that they are obtained by 
the DD method. The matrix Q' contains the sensitivity derivatives of the field variables; the matrix X' contains 


the grid-sensitivity terms (which typically are obtained by differentiating the grid-generation code). The very large 
linear system (Eq. (2.4)) is solved first for Q' so that the SD °F' can be calculated subsequently. 


2.1.2 The Adjoint- Variable (AV) Method. As an alternative to solving Eq. (2.4) for Q', an adjoint-variable 
matrix A is introduced to combine Eqs. (2.3) and (2.4); the matrix A is then specified to ensure that the resulting 
coefficient of Q' vanishes. The AV method becomes 


A p' 


- 9 F y'x 8F 

= 9X X+ 9F + A Ux*' 


9R\ 
9b / 



(2.5) 

( 2 . 6 ) 


The matrix A F' contains the sensitivity derivatives of interest; the superscript A denotes that they are obtained 
by the AV method. However, D F'= A F'= The very large linear system of Eq. (2.6) is first solved for A in 
order that the SD A F' can be calculated subsequently. 

The dimension of b and, thus, the column dimension of Q' is the number of design variables (NDV); the 
dimension of F and, thus, the column dimension of A is the number of output functions (NOF). Therefore, if 
the NDV is greater than the NOF, then the solution of Eq. (2.6) is likely to be computationally less expensive 
than that of Eq. (2.4). (It will definitely be less expensive for a direct-solution procedure; iterative methods are 
normally required, however, because of the extreme size of the coefficient matrix.) 


2.13 The Incremental Iterative Solution Method. As an alternative to pure Newton iteration, typical CFD 
codes employ what is sometimes called quasi-Newton iteration, which is an IIM, to solve the nonlinear flow 
system (Eq. (2.1)); this can be expressed as 

_£51aq" = R" (2.7) 

9Q 

Q1+1 _ QA + A Q». n — 1,2,3, ... (2.8) 
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The left-hand-side (LHS) coefficient matrix operator ^ of Eq. (2.7) is, in many CFD codes, at best only a 
rough approximation to the exact Jacobian matrix operator that is associated with true Newton iteration. Thus, 
Eqs. (2.7) and (2.8) are intended to represent a broad spectrum of implicit and explicit iterative algorithms that 
are common to CFD software. 

Numerous computational difficulties are associated with solving the FO linear sensitivity equations in the 
standard form given by either Eq. (2.4) (the DD method) or Eq. (2.6) (the AV method); these difficulties are 
documented, for example, in Refs. [1], [14]. Previous studies (Refs. [l]-[5]) have shown that these computational 
difficulties can be overcome (at least in part) by iteratively solving these equations in incremental iterative form. 
For the DD method (Eq. (2.4)), the IIM is 

AQ' m = R' m (2.9) 

9 Q 

Q' m+1 = Q' m + AQ' m ; m = 1,2,3,... (2.10) 


where 


Tj/m _ £5;rym , £5:x' + 

R " 9q q + ax x + at 


( 2 . 11 ) 


In Eq. (2.9), the LHS coefficient matrix f|| represents any convergent, computationally convenient approximation 
of the exact Jacobian matrix. In particular, the identical approximate LHS operator and algorithm that are used 
to solve the nonlinear flow equations can also be used to solve the linear sensitivity equations. Comparison of 
Eqs. (2.7) and (2.8) with Eqs. (2.9) and (2.10) reveals that the sensitivity equations are solved by interchanging 
the right-hand side (RHS) of Eq. (2.7) with that of Eq. (2.9) and “freezing” the LHS operator at the steady-state 
value. At convergence, the accuracy of the SD is not compromised if the terms on the RHS of Eq. (2.9) are 
evaluated consistently. The use of the IIM is also applicable in the AV method to solve Eq. (2.6); in this case, 
the LHS operator fg must be transposed. The IIM for this becomes 
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^ AA m = G m (2.12) 

A m+1 = A m + AA m ; m = 1,2,3, ... (2.13) 

where 



and the superscript T indicates a matrix transpose. 

2.2 Applications of ADIFOR 

This section describes different applications of AD to assist in the efficient, accurate calculation of FO SD 
from advanced CFD codes. In particular, the AD precompiler software tool ADIFOR (Automatic Differentiation 
of FORt ran) of Refs. [9] and [15]-[18] is used in this study. The ADIFOR precompiler tool is applied to the 
original FORTRAN program source code from which the SD are to be obtained; the output of this precompiler 
procedure is a new, differentiated source code, which upon compilation and execution will compute (exacdy) 
the numerical value(s) of the derivative(s) of any specified output fiuiction(s) with respect to any specified input 
variable(s). In addition, the new program will perform the function evaluations of the original code. 

2.2.1 Black-Box Applications. Application of ADIFOR to FORTRAN coding of an iterative solution algo- 
rithm (e.g., CFD software) produces a similar iterative algorithm for computing the exact derivatives. However, 
as noted in Refs. [2], [9]-[13], and [19], this latter iterative algorithm obtained from a straightforward “black 
box” (BB) AD application may be neither computationally efficient nor robust; in general, it is not in incremental 
iterative form, even if the original solution algorithm was in that form. The previous BB applications of ADIFOR 
to advanced CFD codes (Refs. [9]-[12]) produced iterative algorithms for SD calculations in which the entire 
flow-solution algorithm is differentiated. 

From the discussion in Refs. [9] and [10], this process whereby the SD are iteratively calculated (following 
the BB use of AD) can be represented conceptually by first combining Eqs. (2.7) and (2.8) (i.e„ the basic CFD 
flow-solution procedure) to yield 

Q n+1 = Q" - P" R“; n = 1,2,3,... (2.15) 

where P n = • Differentiation with respect to b then yields the result 

Q' n+1 = Q'" - P" R.' n - P' n R"; n = 1,2,3,... (2.16) 
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In contrast with Eq. (2.16), a hand-differentiated (HD) implementation of the DD method for FO SD can be 
expressed by combining Eqs. (2.9) and (2.10) (an IIM) to yield 

Q ,m + 1 = Q' m _PR' m ; m = 1,2,3,... (2.17) 

Symbolically, of course, Eqs. (2.16) and (2.17) are equivalent only at convergence of the flow solution, when R n 
vanishes (which also results in the disappearance of P' n R n ), and P n becomes constant at the steady-state value. 
Computationally, however, Eq. (2.17) is potentially much more efficient for several reasons which are listed below. 

1) If the AD-enhanced CFD flow code is executed after the original code has produced a well-converged flow 
solution, then the convergence rate of Eq. (2.16) might be accelerated somewhat. However, differentiation of the 
complete CFD solution algorithm and repeated calculation of its derivatives (represented by P' n in Eq. (2.16)), 
although unwanted and unnecessary, is not avoided. The computationally wasteful, repeated calculation of P ,n 
is likely a significant part of the total work represented by Eq. (2.16). 

2) In Eq. (2.17), the term P is constant. Thus, in principle, only a one-time calculation is required; thereafter, it 
should be stored in memory and reused repeatedly for all iterations. (For 3-D CFD calculations on large grids, the 
computer memory of currently available supercomputers might be too small to store the complete P; obviously, 
in this case P cannot be frozen in memory and reused.) In Eq. (2.16), however, P" (the CFD flow-solution 
algorithm) is updated at each iteration. 

3) The AD-enhanced CFD code will continue to iterate on the solution to the nonlinear flow equations, regardless 
of whether or not they are already well converged. 

4) With the BB application of AD represented by Eq. (2.16), all parts of the term R ,n are forced inside the 
iteration loop and, thus, are calculated at each iteration. However, for the HD IIM represented by Eq. (2.17), most 
of the terms of R ,m can be placed judiciously outside the iteration loop; from Eq. (2.1 1), only the matrix-matrix 
multiplication operation (§§)(Q /m ) must be inside this loop. 

5) The “vectorization” properties of the AD-enhanced CFD code (Eq. (2.16)) for efficient operation on Cray- 
type computers may be severely degraded in comparison with those of the original code. In addition, depending 
on the approach used in the application of the ADIFOR tool and on how many SD are calculated simultaneously, 
the computer memory requirements could become excessive. 

Certain BB applications of AD, discussed subsequently, may enable the complete elimination of the second 
computational difficulty discussed above and greatly limit the impact of the first. Some CFD codes, particularly 
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those of the 2-D implicit type, are equipped with an optional computational-work (CW) saving strategy; it is 
known as the “frozen Jacobians” (FJ) option. This scheme takes advantage of the fact that as the quasi-Newton 
flow-solution method of Eqs. (2.7) and (2.8) converges, the LHS operator of Eq. (2.7) becomes approximately 
constant The FJ option provides the capability of freezing (not updating) these terms (represented by P n in Eq. 
(2.15)) for a specified number of iterations; the result is typically a large savings in CW per iteration. 

For an AD-enhanced CFD code with the FJ option (henceforth known as the BBFJ method), the potential for 
CW savings is very large (i.e., proportionally far greater than for the original CFD code because the unnecessary 
repeated update of P n and of the unwanted P' n can be avoided in Eq. (2.16)). The AD-enhanced CFD code 
is simply started with a well-converged flow solution, and the FJ option is set to activate for all iterations (after 
the first iteration). 

With the BBFJ strategy, clearly the first computational inefficiency (discussed previously) is reduced 
significantly, but not eliminated. It is suggested that further improvements to the BBFJ method might be made 
by editing out these terms (i.e., P /n R n of Eq. (2.16)) from the AD-enhanced code. If successful, this process 
could also have significant collateral benefits with respect to reduced computer memory requirements and restored 
vectorization. These improvements are aimed at making the BB method (Eq. (2.16)) more like the HD DD IIM 
(Eq. (2.17)). 

2 . 2 3 . Incremental Iterative Applications. Earlier studies (Refs. [2] and [9]-[ll]) have proposed that many 
of the previously discussed computational inefficiencies associated with the BB application of ADIFOR to CFD 
codes (Eq. (2.16)) can be overcome (at least in part) by a more judicious application of ADIFOR. The goal of 
this approach is that the resulting SD calculations are (more nearly) in the IIM form of a HD application of the 
DD method (Eq. (2.17)). Specifically, in the present study, ADIFOR is applied to differentiate only the RHS 
of Eq. (2.7), which is the residual R of the nonlinear flow equations (Eq. (2.1)). These differentiated terms are 
assembled on the RHS of the IIM of Eq. (2.9); thus, the resulting scheme is essentially that expressed by the 
efficient Eq. (2.17) which is an IIM. The construction of the required derivatives is now via AD rather than HD. 
This scheme is henceforth known as the ADII method; it should effectively combine an existing, highly efficient, 
iterative solution algorithm with a fast, accurate, reliable procedure for constructing all required terms. 

The ADII strategy will bypass the most obvious computational inefficiencies of the BB strategy. For example, 
the unnecessary construction and repeated evaluation of the term P /n in Eq. (2.16) is completely avoided, and the 
inverse approximate operator P is not updated at each iteration (at least in principle). Evaluation of all derivative 
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terms except |gQ' m can be placed outside the iteration loop; however, note that the HD evaluation of fgQ ,m 
would be less costly (recall that only the explicit matrix-matrix multiplication is then required at each iteration). 
Although the repeated calculation of R n is not avoided, the repeated full iteration on the nonlinear flow equations 
does not continue with the ADH scheme. Despite these improvements, some important issues remain in regard 
to the vectorization properties and computer memory requirements associated with the ADH procedure. These 
issues are partially addressed herein; see Ref. [20] for expanded discussion. 

A more detailed discussion is provided in Ref. [20] of how ADIFOR is applied to implement the ADII 
scheme; only a few key highlights are mentioned here. Most importantly, this scheme must be assembled with 
great care to ensure that contributions to the SD from the boundary conditions are taken into account fully; 
failure to do this will result in severe errors in the calculated SD (Ref. [21]). In the present study, this involves 
separate applications of AD (i.e., applications to a master boundary-condition subroutine and applications to a 
master interior-cell-residual subroutine). Furthermore, these applications of AD are further subdivided to ensure 
the terms that must be calculated inside the iteration loop (i.e., can be separated from the remaining 

terms that should be placed outside the loop. Finally, the AD versions of these subroutines are then carefully 
interwoven to function as the ADD scheme*. 

One important and useful feature of the ADIFOR system for AD is that terms of the type f^Q' or f^X' 
(recall Eqs. (2.4) and (2.11), for example) are calculated without the explicit calculation of the very large Jacobian 
matrices or |§, respectively, and without explicit postmultiplication by the matrices Q' or X', respectively. 
Of course, the AD-enhanced code, which can evaluate these complete expressions, will require increased memory 
over that of the original code. However, this increase is approximately equal only to the memory of the original 
code times the column dimension of Q' and X'. For the present application, this is NDV, which is the dimension of 
b (or the dimension of that fraction of b for which SD are to be concurrently calculated in the ADII method). The 
final result is an extremely fortuitous conservation of computer memory; without this conservation of memory, 
given the overwhelming size of and §§, the application of ADIFOR to advanced CFD codes would be 
infeasible. 

In contrast with the preceding discussion, expressions of the form (§§) X (recall Eqs. (2.6) and (2.14) of 
the AV method) cannot currently be evaluated via applications of ADIFOR without the explicit calculation of 
the very large transposed Jacobian matrix (§§) and the postmultiplication of it by the matrix A. For modern 

* This demonstration of ADII is for a single grid code; the ADII scheme is currently under development for a multigrid, multiblock code. 
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CFD codes, this step is completely infeasible. Unfortunately, ADIFOR cannot be considered to assist in the 
construction of these terms in the IIM for the AV scheme (Eq. (2.14)); nevertheless, it can be used to construct 
the remaining terms (Eq. (2.5)). 


2.2.3 Turbulence-Modeling Applications. The presence of turbulence modeling presents a challenge in the 
calculation of aerodynamic SD. Practically speaking, the task of differentiating the turbulence-modeling terms 
to include their influence in the Jacobian matrices and other terms of the DD and AV methods is too complex 
to do by hand. Symbolic manipulators could be used to differentiate the algebraic equations involved, with 
program-flow control by macros, and then to generate SD code; however, ADIFOR has the advantage of being 
able to directly work with the existing FORTRAN source code, with automatic program-flow control and global 
dependency checking. 

In the earlier study of Ref. [1], a HD IIM version of the DD and AV schemes was created to compliment a 2-D 
TLNS CFD code. (Henceforth, these two HD schemes are referred to as the DDH and AVII methods, respectively.) 
These two schemes were shown to generate very accurate SD for constant-viscosity laminar flow but produced 
significantly erroneous SD for turbulent flow; this discrepancy resulted because the turbulent viscosity terms (from 
the Baldwin-Lomax algebraic model (Ref. [22]) were not differentiated by hand because of their complexity. In 
the present study, ADIFOR is applied to correct this deficiency. That is, ADIFOR is applied to differentiate 
the turbulence-modeling terms only, and the results are incorporated as a correction in the HD SD code. This 
correction is applied to the DDII scheme only, which results in a method known henceforth as DDIITC. This 
correction could not be added in full to the AVII scheme, for reasons discussed subsequently; therefore, no 
AVHTC strategy currently exists. 


Conceptually, the AD correction for the viscosity terms is added to the DDII scheme as 

9F _ /9F\ 9F 9V 
9Q ~ V9Qyt + 9V 9Q 
9F _ /9F\ 9F 9V 

ax ~ {dxl + avax 

9F _ / 9F\ 9F 9 V 
9b “ V9b/l + 9V 9b 


(2.18) 

(2.19) 

( 2 . 20 ) 


where V is a vector of viscosity terms (including the turbulent viscosity). The subscript V in the above indicates 


differentiation within the term with V held constant; thus, terms with this subscript represent the original terms 


of the uncorrected HD code. Substitution of the above into Eq. (2.3) of the DD method results in 
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where 


v' = 


dV 

db 




Similar manipulations applied to Eq. (2.4) of the DD method yield 


( 2 . 22 ) 




(2.23) 


Clearly, the viscosity-derivative correction terms to be added in Eqs. (2.21) and (2.23) are |^V' and |^V', 
respectively. In the present study, fy- and are constructed by hand, and V' is constructed via ADIFOR. The 
three terms of V' on the RHS of Eq. (2.22) are constructed with separate applications of ADIFOR so that the 
terms f^Q' can be separated from the others and placed inside the iteration loop. The terms fgQ' and f^X' are 
assembled without the explicit computation of the Jacobian matrices fg and f^. Thus, these terms are evaluated 
without an excessive expansion of the computer memory, as discussed previously in regard to the AD-assisted 
evaluation of the terms §f|Q' and f^X'. 

The DM for solving Eq. (2.23) becomes 



Q' m+1 = Q' m + AQ' m ; m = 1,2,3,... 


(2.24) 

(2.25) 


where 




V' m = — — Q' m 


av , av 
+ ax x + ab 


(2.26) 

(2.27) 


In the turbulent sample problem of this study, the term |^Q /m was a computationally expensive addition 
to the iteration loop, even after the ADIFOR-generated code was extensively massaged to restore vectorization 
and other features related to efficiency. Initially, the computational cost per iteration was about 3.62 times more 
costly per iteration when the correction was switched on, although the overall rate of convergence was not affected 
greatly. A CW saving strategy was proposed and tested, where the term §gQ' m of Eq. (2.27) was frozen (not 
updated) inside the iteration loop for a specified number of iterations. For 10 frozen iterations prior to each update 
of this term, the overall increase in average CW per iteration due to this turbulence-modeling correction was about 
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26.6 percent (compared with the CW per iteration with the correction switched off) with no major impact on the 
rate of total error reduction. Furthermore, the correction had little impact on the computer storage requirements. 
The unsuccessful attempt to correct the AVII scheme for turbulent flow yielded the following formulation: 


'-(SMSMCSMSl] 

♦(5 ♦**)(**♦£) 


(2.28) 


-di' 


AA m = G m = G^ + G)J? C 
A m+1 =A m + AA m ; m = 1,2,3, ... 


(2.29) 

(2.30) 


where 

g " - (Hl A " + (SX 

All parts of the turbulence correction can be constructed easily via ADIFOR, except G!p c (Eq. (2.32)). The 
present version of ADIFOR can construct this term only by explicitly computing the large matrix and 

postmultiplying it by the terms shown; as discussed previously for the matrix-matrix product (§§)a, this 
procedure is not feasible for modern CFD codes. 


2.2.4 Vectorization and Memory Considerations. Prior to the compilation and execution of any AD- 
enhanced FORTRAN source code, a parameter g$p$ is specified within the code. For each execution of the 
code, this parameter determines the number of independent (design) variables with respect to which derivatives 
are concurrently computed. Thus, the user has the following options: 

1) Compute all required derivatives by executing the AD-enhanced code once for each independent variable 
(i.e., NDV code executions with g$p$ = 1). 

2) Compute all required derivatives by executing the AD-enhanced code only once (i.e., one code execution, 
with g$p$ = NDV). 

3) Set g$p$ such that 1 < g$p$ < NDV; this requires multiple executions (less than NDV) of the AD-enhanced 
code, where subgroups of g$p$ derivatives are concurrently computed for each code execution. 
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The specified value of g$p$ has a significant impact on computational requirements in several critical ways. 
With respect to memory, for example, recall that the memory increase of the AD-enhanced code is approximately 
equal to g$p$ times the memory of the original code; thus, if this parameter is too large, the memory requirements 
of the code could be excessive. 

The AD-enhanced code retains all do-loops and function evaluations of the original code. Within each original 
do-loop is inserted one or more new innermost do-loops; the length of each new do-loop is g$p$ (e.g., DO 10 
1=1, g$p$). Inside these new loops, derivative calculations are made. The presence of these new innermost 
do-loops has a profound impact (frequently negative) on the vectorization characteristics for performance on 
Cray-type computers: 

1) The do-loops of the original code, which previously vectorized, will no longer be vectorized in the 
AD-enhanced version. An exception to this is when g$p$ < 5; the “aggressive” Cray compiler option will 
automatically “unwind” the new innermost loops and may restore the vectorization of the original loops, complete 
with the derivative calculations. 

2) For g$p$ > 6, vectorization of the original loops is not recovered, but with the aggressive compile option, 
the new innermost loops are vectorized. Nevertheless, overall code performance remains poor on Cray computers 
unless g$p$ is large enough that the vector lengths become sufficiently long for efficient execution on these 
machines. At the same time, however, for large g$p$ the computer memory requirements of the AD-enhanced 
CFD software can become excessively large. 

Apart from the vectorization considerations discussed above, the number of arithmetic operations per 
concurrently computed derivative is always decreased as g$p$ increases. This happens because, for each execution 
of an AD-enhanced code, part of the derivative calculations occur outside of the innermost loops, and the results are 
reused for all derivative calculations within the innermost loops. Furthermore, the complete function evaluations 
of the original code are performed. 

The sample problems illustrate the consequences discussed previously; the results from these sample problems 
are to be given. For example, (except when g$p$ is large) g$p$ = 5 produces the highest computational efficiency 
per design variable, and this efficiency is progressively reduced as g$p$ is reduced to 1. A particularly inefficient 
case is that of g$p$ = 6 (thereafter efficiency gradually increases as g$p$ increases). In the case with NDV = 6, 
rather than perform one code execution with g$p$ = 6, two code executions, each with g$p$ < 6 (e.g., the first 
execution with g$p$ = 5 and the second execution with g$p$ =1), were significantly more efficient. 
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23 Computational Results: FO SD 


1\vo sample problems are considered here; they are identical in every way to those studied previously in 
Ref. [1], where a more complete description is given. The first example is low-Reynolds-number (Re = 5x 10 3 ) 
subsonic (Moo = 0.6), constant-viscosity laminar flow over an isolated NACA 1406 airfoil at an angle-of-attack, 
a = 1.0°. The second example is similar, except the flow is a high-Reynolds-number (Re = 5x 10 6 ) transonic 
(Moo = 0.8) turbulent flow. Flow calculations are made on a C-mesh with dimension 257x65 (circumferential 
x normal direction); the clustering of points near the airfoil’s surface was tighter for the high-Reynolds-number 
example. Grid sensitivity derivatives were produced with a unique scheme that was first reported in Ref. [14] and 
was subsequently applied to these sample problems in Ref. [1]. 

The CFD code applied here (and in Ref. [1]) solves the 2-D TLNS equations with an upwind, cell-centered 
finite-volume formulation with a higher order accurate evaluation of all fluxes and the algebraic turbulence 
modeling of Baldwin and Lomax (Ref. [22]). The code employs an implicit, spatially split approximate- 
factorization flow-solution algorithm. Also available is the FJ option (discussed previously), where the entire 
implicit operator (i.e., the complete set of LU-factored block-tridiagonal coefficient matrices) is stored in memory 
and repeatedly reused (not updated) for a specified number of iterations. 

As in Ref. [1], the FO SD of three aerodynamic output functions Cl, Cd, and Cm (the coefficients of 
lift, drag, and pitching-moment, respectively) are calculated with respect to three geometric shape variables T, 
C, and L (maximum thickness, maximum camber, and location of maximum camber, respectively) and with 
respect to three flow variables a , M TO , and Re (each defined previously). Therefore, F = (Cl, Cd,Cm) and 
b = (T, C, L, a , Moo, Re) T . The SD are computed with a wide variety of different methods, and the results are 
compared on the basis of accuracy and computational time and memory. 


23.1 Accuracy Comparisons. The FO SD are calculated for both sample problems with the methods CD, 
DD1I, AVII, ADII, and BB; in addition, the DDIITC scheme is applied only to the turbulent example (because 
it is unnecessary for the laminar case). The application of the CD, DDR, and AVII schemes to these problems 
repeats the work of Ref. [1]; the manner in which these schemes are applied is discussed in depth in Ref. [1], Of 
course, the DDIITC, ADII, and BB schemes are the methods for which derivatives are calculated via applications 
of ADIFOR, either in part or in total (depending on the scheme). 
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The SD that were calculated are presented in Table (A.l) (in Appendix A) for the laminar example. The 
actual numerical values of the SD are given for the CD method. For the other methods, SD ratios are given (i.e., 
each SD has been normalized by the respective SD calculated via the CD scheme). Table (A.l) clearly shows 
excellent agreement among all these methods, as expected. 

The SD for the turbulent example are presented in Table (A.2). The actual SD are shown for the CD 
scheme; the remaining cases are shown as SD ratios. As expected, this table shows excellent agreement within 
the DDHTC, ADII, BB, and CD methods. The results for the DDII and AVII methods do not agree well (for 
some SD) with the other schemes because of the turbulence-modeling terms, as discussed previously herein and 
in Ref. [1], The erroneous results for the DDII and AVII schemes agree extremely well with each other, however, 
because the two are algebraically equivalent; for this reason, the results for these two methods are shown as a 
single result in Table (A.2). 

The good accuracy and agreement in the preceding results is due in part to the very tight convergence 
tolerances that were enforced on all calculations. The average total error was reduced to machine zero (a relative 
reduction of approximately 12 orders-of-magnitude (OM)) in the initial flow solution and in all twelve flow 
solutions that were required for the CD method (i.e., two solutions per design variable). Very small forward and 
backward perturbations Abj= ±5.013-6 x bj are made to each design variable to ensure good accuracy with the 
CD method. For each linear system that was solved to compute these SD, the error was reduced at least eight OM; 
in each sample problem, this involved NDV = 6 solutions for the DDII, DDHTC, ADH, and BB methods, and 
NOF = 3 solutions for the AVII scheme. Of course, these tight convergence tolerances are far more restrictive than 
would be required in ordinary engineering practice and greatly increase the computational cost of the calculations. 

2.3.2 Computational Time and Memory Comparisons. In this section, some of the methods discussed 
in the previous section are further subdivided; these subdivisions have little or no impact on the SD that are 
calculated but can have a significant impact on the total computational efficiency of the method. 

The CD method is subdivided into two methods, depending on whether or not the FJ option is activated; 
when active, 10 iterations of Eqs. (2.7) and (2.8) are specified prior to each update of the LHS operator. The 
methods without and with the FJ option are referred to as the CD and CDFJ methods, respectively. 

Similarly, the BB method is subdivided into two methods, depending on whether the FJ option is activated 
(the BBFJ scheme) or not (the simple BB scheme). These two BB methods are further subdivided into additional 
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methods, depending on how the parameter g$p$ is applied. Two options were tested: first, all SD were calculated 
with a single execution of the AD-enhanced code (i.e., this implies g$p$ = NDV = 6); and second, two executions 
of the code were made — the first with g$p$ = 5 and the second with g$p$ = 1. From these subdivisions, the 
methods are BB(6), BB(5+1), BBFJ(6), and BBFJ(5+1). Execution of the AD-enhanced code is started for all of 
these BB-type methods after the original code has produced the fully converged flow solution. (Recall that the 
fully converged flow solution is required at start-up for the latter two methods.) 

The ADII scheme also is subdivided into two methods, depending on the application of g$p$ to yield the 
ADII(6) and ADII(5+1) methods. For turbulent flow, the DDHTC method is further subdivided according to 
whether or not the option to freeze the turbulence-correction terms is activated. The terminology DDIITCFR 
indicates that this option is activated. (Ten frozen iterations are specified for each iteration that updates these 
terms.) The notation DDIITC indicates the option is not activated. 

Comparisons of the time (total CPU time) are shown in Table (A.3) for both the laminar and turbulent sample 
problems; all calculations are performed on a Cray-YMP computer; results are given in secs. All reported timings 
do not include the cost of the initial flow solution. Note that a superscript * in the tables indicates estimated 
results, based on results in Table (A.4) which shows a comparison of the computational times (the CD and CDFJ 
methods are excluded) in CPU time per iteration per linear system solved. The results from Table (A.3) for 
the methods ADII(5+1), BB(5+1), and BBFJ(5+1) have been separated in Table (A.4) to compare the individual 
effect of g$p$ = 5 and g$p$ = 1. 

The total memory requirements for the different methods are compared in Table (A. 5), where the results 
are given in Mega- words (Mw); no difference occurs in the results for the laminar and turbulent examples in 
this table. The computer memory requirements for the ADII scheme are less than that for the BB approach, 
particularly if the original flow code is of the type that uses a large amount of memory for storage of the terms 
of the LHS operator; recall that these terms are not differentiated with the ADII scheme, which results in a 
significant conservation of computer memory. 

2.4 Conclusions: FO SD 

Conclusions based on the calculations for FO SD are enumerated subsequently. 

1) The HD IIM schemes, although presently the most efficient, are very difficult and time consuming to construct 
accurately, even for relatively simple CFD codes. For more complex codes with features such as turbulence 
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modeling, etc., this approach is not feasible. ADIFOR is a reliable tool for the quick construction of accurate 
source code to evaluate all or parts of the SD from complex CFD codes, but straightforward BB application of 
ADIFOR to CFD codes can be slower than even the CD method and require substantially more memory than 
the original code. 

2) ADIFOR can be used successfully to create corrections to HD SD codes, where only relatively small, previ- 
ously undifferentiated parts of the original flow code (such as turbulence-modeling subroutines) are differentiated 
via ADIFOR. These corrections can be very costly with respect to the efficiency of the HD SD code, as seen 
in Tables (A.3) and (A.4) for the turbulent-flow problem; the cost in computer memory is negligible, as shown 
in Table (A.5). 

3) For all applicable methods, the computational penalty associated with the turbulence-modeling terms (all 
constructed via ADIFOR) is significant and disproportionately high. Detailed comparisons of the laminar and 
turbulent timings (for a given method) shown in Table (A.4) reveal this penalty. The disproportionate cost is 
highest for the most efficient method and also for the application where g$p$ = 1. The AD-correction for the 
turbulence-modeling terms in the DDIITC method represents an inefficient (g$p$ = 1) application of ADIFOR; 
the inordinately large computational cost of this turbulence correction is thereby explained. 

4) The ADH scheme is not as efficient, computationally as the HD IIM schemes but is more efficient than all 
other methods tested. For example, depending on the particular sample problem and application of g$p$, the 
ADII scheme produces computational improvements by a factor which varies from approximately 6 to 15 when 
comparisons are made with the simple BB scheme; similar comparisons to the BBFJ scheme yield factors which 
vary from approximately 1.4 to 1.7. The most efficient ADII results were approximately two to four times more 
efficient than the results from the CDFJ scheme — even more efficient compared with the CD method. 

5) The ADII scheme is not as easy to implement as the BB methods. For example, particular care must be 
talfpn to ensure that the contributions from the boundary conditions are properly taken into account. However, 
when compared with the HD approach, the ADH scheme can be implemented easily with very accurate results, 
even for very advanced CFD codes. For example, the time required to develop the source code for some of these 
different methods is estimated; HD IIM (DDII, AVII, etc.)— six man-months to two man-years, or even longer, 
depending on the complexity of the flow code; ADII — about one man-week; BB — about one man-day. 

6) The BBFJ strategy is no more difficult to implement than the simple BB approach, if the original flow code 
is equipped with the FJ option. A very large increase is noted in the computational efficiency (compared with the 
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simple BB strategy) when this option is activated; this improvement was by an impressive factor which varied 
from approximately four to nine, depending on the particular example problem and application of g$p$. The 
BBFJ strategy is not as efficient as the ADII scheme with respect to computational time or memory. However, 
with relatively minor code modifications, the method possibly could be made to function with nearly the efficiency 
of the ADII scheme. Therefore, when the FJ option is available, the BBFJ scheme is the simplest method to 
implement which also maintains reasonable efficiency (compared with the CDFJ scheme in particular); otherwise, 
the necessary extra effort should be invested to implement the ADII scheme. 

7) For the BB, BBFJ, and ADII schemes, the computational cost in terms of CPU time and computer storage is 
very sensitive to the value selected for the parameter g$p$; the effect can vary significantly for different machines. 
The significance of g$p$ is critical for a large NDV; on Cray computers with the aggressive compile option, the 
choice of g$p$ = 5 seems to provide the most efficiency in terms of CPU time per iteration per linear system 
solved, with a manageable increase in memory. 

8) Currently, ADIFOR cannot be applied to construct the AV method in total; this limitation is a serious 
consideration because great potential exists for efficiency with the AV scheme when NDV is much larger than 
NOF. 


3.0 Second-Order Sensitivity Derivatives (SO SD) 

3.1 Basic Equations and Incremental Iterative Forms 

A brief derivation is presented of the basic equations of SO discrete aerodynamic sensitivity analysis. The 
result is four methods, denoted as (1) DD.DD, (2) AVDD, (3) DDAV, and (4) AVAV; this notation roughly 
parallels the derivation and description given in Ref. [23] for SO shape sensitivity analysis applied to linear 
heat-conduction problems. In addition, the incremental iterative forms are given for solving the additional large 
linear systems that result from Methods (1) and (2). 

For convenience and subsequent notational clarity, the key equations for the FO derivatives are repeated, 
where only the terms for the i-th aerodynamic output function (F,) and for the j-th design variable (bj) are given 
here. Recall the FO DD approach is 


d f , - 4F S _ 3Fi 9Fi . 9^ 

,;J - dbj aQ^ + ax A ' + abj 

(3.1) 

dR_ 8R aR x , aR_o 

^“dbj 6Q Q) + aX Xj + abj 

(3.2) 
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where Qj = Xj = ji|X. The FO AV approach is obtained by introduction of the adjoint-variable vector Ai 
into the preceding equations to eliminate Q'; the result is 


fi;j = 



(3.3) 

(3.4) 


The following differential operator is defined for use in the derivations that follow. 

D( ) _ 8( ) dQ a( ) dx 9( ) 

Dbfc 9Q dbjj 9X dbjj 


Therefore, for example, comparison of this operator with Eqs. (3.1) and (3.2) yields, respectively. 


D(Fj) _ D , . 
Dbj * iJ ’ 


P(B-) _ |j> 
Dbj ’ 


(3.6) 


3.1.1 Method 1: DD.DD. An inspection of Eqs. (3.1) and (3.2) reveals explicit dependencies in each on Q' 
and Xj, in addition to Q, X, and b; this dependency is expressed as 


TT 7 = D Fj ; j = D F/.j (Qj (b) , Xj f (b) ; Q(b), X(b), b) 


db 


dR 

db, 


= Rj = Rj (Q](b)t Xj(b); Q(b),X(b),b) 


Differentiation of Eqs. (3.7) and (3.8) with respect to the k-th design variable b k yields, respectively, 

d 2 Fj _ 9Fi 9F; , D ( PF i;i) 

db k dbj 9Q Wj,k 9X J ' k Db k 

d 2 R _ 9R 9B. „ D ( R l) _ 

k. jk. Air j> k + Db k 


(3.7) 

(3.8) 


(3.9) 

(3.10) 


db k dbj 9Q 

where Qj' k = dbk J b . and X j k = dbkdbj • 

In the preceding differentiation, the chain rule is applied term by term to Eqs. (3.7) and (3.8). Equations 
(3.1) and (3.2) are used to produce the simplifications that result in the first two terms of Eqs. (3.9) and (3.10), 
respectively; each of the third terms is very complex and has been simplified as a single term with the special 
differential operator defined by Eq. (3.5). A more complete expansion of these terms is provided in Ref. [20]. 
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Clearly, Xj' k is the SO grid-sensitivity term, which is obtained in general (for geometric design variables 
only) by twice differentiating the mesh-generation code. These SO grid-sensitivity terms vanish in the example 
problems of the present study because of the linear nature of the particular remesh/grid-sensitivity scheme used. 
(See Refs. [1] and [14].) 

In all subsequent discussions, it is assumed that the complete Hessian matrix ^ is desired for each output 
function F;. The DDDD method requires a priori knowledge of the complete Q' matrix, which is obtained via 
NDV solutions of the FO DD system (Eq. (3.2)). Thereafter, a maximum of (NDV) 2 solutions of the SO system 
(Eq. (3.10)), is required to determine all Qj' k ; however, if the identity Q" k =Q kj is exploited computationally, 
then the minimum number of SO solutions is [(NDV) 2 + NDV]/2. Thus, the DD.DD method requires a minimum 
total number of NDV + [(NDV) 2 + NDVJ/2 solutions of very large linear systems. 

The coefficient matrices of the FO DD system and the SO DDJDD system are identical; thus, when these 
systems are cast in incremental iterative form, both could be solved with the identical approximate LHS operator 
and algorithm that is also used to solve the nonlinear flow equations. The IIM for solving the SO Eq. (3.10) is 


AQ\ m = 
9Q W)k 


ft m u// in 

R j>k 


(3.11) 


Q" m+1 


Q h m i An" re 

i,k +^Wj,k 


m = 1,2,3,... 


(3.12) 


where 


n" m _ ^ rV* m 

R l> - 


* £5x" * 

+ ax j ' k+ Db k 


(3.13) 


3.1.2 Method 2: AV.DD. An inspection of Eq. (3.3) reveals explicit dependencies on Aj and Xj; in Eq. 
(3.4) explicit dependence on A; occurs but not on Xj. Both equations depend explicitly on Q, X, and b. The 
complete dependencies are expressed as 

= A F i ' ;i = A F i ' ;j (A; (b), Xj(b); Q(b), X(b), b) (3.14) 

abj 

Gi= Gi(Aj(b);Q(b),X(b), b) (3.15) 


Differentiation of Eqs. (3.14) and (3.15) with respect to bk yields, respectively, 

d 2 Fj 


db k db; 


i . (W .aT^x" 

b • = (ax x ‘ + ^ j Ai;k + lax + ax J x i> 


+ 


p ( A ni) 

Db k 


(3.16) 
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dGj 

db„ 


( 3 . 17 ) 



D(Gj) 

Db k 


= 0 


where A?. k = 

In the preceding differentiation, the chain rule is applied to Eqs. (3.14) and (3.15). Equations (3.3) and (3.4) 
are used to produce the simplifications that result in the first two terms of Eq. (3.16) and the first term of Eq. 
(3.17); each of the last terms is complex and has been simplified as a single term with the operator defined by 
Eq. (3.5). (See Ref. [20] for a more complete expansion of these terms.) 

The AV.DD method requires a priori knowledge of the complete Q' and A matrices, which are obtained via 
NDV solutions of the FO DD system (Eq. (3.2)) and NOF solutions of the AV system (Eq. (3.4)), respectively. 
Thereafter, NDV x NOF solutions of the SO system (Eq. (3.17)) are required to determine all A;. k ; the AV.DD 
method thus requires a total of NDV + NOF + (NDV x NOF) solutions of very large linear systems. 

The coefficient matrices of the FO AV system and the SO AVDD system are identical; thus, when these 
systems are cast in incremental iterative form, both could be solved using the identical approximate LHS operator 
and algorithm (the transpose of that which is also used to solve the nonlinear flow equations). The IIM for 
solving the SO Eq. (3.17) is 


/ \ T 

/ ^ | a a * m — r* * m 

-[ sq ) AA i*- G i;« 






( 3 . 18 ) 


where 



( 3 . 20 ) 


3.13 Method 3: DD.AV. This method is derived by introducing an arbitrary adjoint-variable vector into the 
DD£>D method to combine Eqs. (3.9) and (3.10); the adjoint-variable vector is specified so that the resulting 
coefficient of Qj' k vanishes. The resulting DDAV method is 


d 2 Fj 

db k dbj 


/ dFj -t 

Ux +Ai axj x i' k 

A r 0 


+aT- 


Db k 


P ( DF . ? ;i) 

Db k 


( 3 . 21 ) 
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The adjoint- variable vector A; in the derivation of Eq. (3.21) is identical to that of the FO AV method and is 
found by solving the FO Eq. (3.4). 

The DD.AV method requires no simultaneous solutions of large systems of linear equations that involve SO 
terms (in contrast with the previous two schemes). An inspection of Eq. (3.21) reveals that the method requires 
knowledge of the complete Q' and A matrices; they are obtained via NDV solutions of the FO DD system 
(Eq. (3.2)) and NOF solutions of the FO AV system (Eq. (3.4)). respectively. Thus, the DD.AV method requires 
a total of NDV + NOF solutions of large simultaneous systems of linear equations; only the systems for FO 
derivatives are solved. 


3.1.4 Method 4: AV.AV. This method is derived by introducing a new, arbitrary adjoint-variable vector into 
the AVX)D method to combine Eqs. (3.16) and (3.17); this adjoint-variable vector is specified so that the resulting 
coefficient of A[ ;k vanishes. The resulting AVAV method is 


d 2 Fj 

db k dbj 


p ( A ni) 

Db k 


+ (Q;') T 


D(Gi) 

Db k 


(3.22) 


The adjoint-variable vector (which results in the disappearance of terms that involve A[ ;k ) is seen in the derivation 
of Eq. (3.22) to be Q-, which is found by solving the FO DD (Eq. (3.2)); detailed demonstration of this result 


is given in Ref. [20]. 


No computational advantage is associated with the AV.AV method over the DDAV method; in Ref. [20], 
the two schemes are shown to be term-by-term equivalent. That is 


(Q i> -5ST 1 


D ( R 0 

Db k 


D ( A Fj ; j) _ D( D F( ;j ) 
Db k Db k 


(3.23) 

(3-24) 


The AV.AV method requires knowledge of the complete Q' and A matrices, which are obtained from solving 
Eqs. (3.2) and (3.4), respectively. Thus, the AVAV method requires a total of NDV + NOF solutions of large 
systems of linear equations; only the systems for FO derivatives are solved. 


3.1.5 Discussion. An analysis was made to determine which of the preceding four methods would be potentially 
the least costly computationally by considering the of total number of large simultaneous linear systems that must 
be solved to calculate the complete Hessian matrix for all Fj. The conclusions of this study are: 

1) Methods (3) and (4) are computationally equivalent and are henceforth known as Method (3/4). 
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2) Method (3/4) is unconditionally less costly than Method (2); therefore, either Method (1) or Method (3/4) 
should be selected, depending on conclusion 3, which is given subsequently. 

3) If the equality Qjj k =Q k j is fully exploited, then Method (1) is less costly than Method (3/4) when NDV 
x (NDV + 1) is less than 2 x NOF. If the same equality is not exploited, then Method (1) is less costly than 
Method (3/4) when (NDV) 2 is less than NOF. 

3.2 Applications of ADIFOR 

Ih this section, various procedures are outlined whereby AD might be applied effectively to assist in computing 
the SO aerodynamic SD. Many of the terms in each of the preceding four methods are exceedingly complex; in 
particular, this applies to the large, complex groups of terms that are symbolized compactly with the operator 
Practically speaking, differentiation and coding by hand to construct these terms is impossible, even for the less 
complicated CFD codes (e.g., 2-D Euler codes). Simply stated, without AD the equations for SO aerodynamic 
sensitivity analysis cannot be constructed. 

3.2.1 Noniterative Applications for Method (3/4) (DD.AV / AV.AV). Fortunately, ADIFOR is ideally 
suited for the quick and reliable generation of source code to accurately evaluate the required SO terms without 
an excessive expansion of the computer memory. For example, the ADIFOR-assisted construction of the SO 
Method (3) (DD.AV, Eq. (3.21)) involves creation of a source code that evaluates » . The source code can be 
created easily by the straightforward application of ADIFOR to an existing subroutine that evaluates Rj; this AD 
application is completely analogous to the AD-assisted creation of source code that evaluates which is Rj, 
from an existing subroutine that evaluates R. (Recall that the creation of source code via ADIFOR to evaluate 
Rj was an important requirement in the success of the FO ADR scheme.) Similar remarks, of course, apply in 
the AD-assisted creation of for use in Eq. (3.21). Therefore, the fact that the a priori calculation of the 

FO Q' matrix is required by all four SO methods (including Method (3)) is somewhat fortuitous because this 
ensures that the FO source code to evaluate Rj (and d F-.j) will be available for further use in the AD-assisted 
creation of the SO terms. 

Except for the terms Xjj k and Aj, the remaining terms in the SO Eq. (3.21) are taken as is from the same 
FO equations that are used to calculate all Qj. The term Xj' k would be produced by twice differentiating the 
grid-generation code, and all Aj must be obtained by solving first the FO AV equations. After the required FO 
equations for the required Q' and A are solved, the SO Method (3) becomes a noniterative, computationally 
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efficient scheme for computing SO SD, where all required SO terms are easily constructed via straightforward 
applications of ADIFOR to key parts of the existing FO source code. The discussion of SO Method (3) in the 
preceding two paragraphs is easily extended to the SO Method (4) because, as noted previously, the two methods 
(Method (3/4)) are equivalent. 

The requirement to first solve the FO AV equations is an important consideration for each SO method 
except Method (1) (DD.DD). Unique computational difficulties exist (discussed earlier in greater detail) which 
are associated with the use of ADIFOR to construct key parts of the AV methods; thus, these difficulties are 
transmitted to these SO methods. This concern speaks to SO Method (3/4) directly, because it was concluded 
earlier that this method is unconditionally more efficient than SO Method (2) (AV ,DD). 

For the particular combination of NDV = 6 with NOF = 3 (which is applicable to the sample problems 
of this study), the SO Method (3/4) would be about 3 times less costly than the remaining choice Method (1); 
this projection is based on the previous discussion, which considers the comparative total number of large linear 
systems that must be iteratively solved (27 large system solutions for Method (1), compared to 9 for Method 
(3/4)). The actual implementation of the SO Method (3/4) is not included in the present study but is a topic 
of ongoing work. 


3.2.2 Black-Box and Incremental Iterative Applications for Method (1) (DDJ>D) and Method (2) 

(AV.DD). The SO Method (1) has been projected to be less costly than Method (3/4) for certain combinations of 
NDV and NOF (which have been specified). An advantage unique to the SO Method (1) is that no adjoint-variable 
equations are to be solved; thus, in principle, the entire scheme can be constructed from start to finish via the 
present version of ADIFOR. However, an important disadvantage also exists: in contrast with the SO Method 
(3/4), large linear systems must be iteratively solved for the SO Q" terms. Most of the remainder of this section 
focuses on ADIFOR-assisted BB and/or UM implementations of this SO method (DDJDD). The discussion and 
resulting methods are analogous to that seen earlier for ADIFOR-assisted implementations of the FO DD scheme. 

In principle, SO aerodynamic SD can be obtained by a BB application of the present version of the AD 
tool (ADIFOR 1.0) to an AD-enhanced version of the original flow code (called the BB.BB method). However, 
the BB.BB scheme can also be obtained by applying the future new version of the AD tool (ADIFOR 2X) to 
the original flow code; this advanced version of ADIFOR is presently being developed with the new, optional 


26 



capability of providing SO derivatives* . For advanced CFD codes, the BB.BB method may be exceedingly 
inefficient computationally, for reasons that are discussed subsequently; in addition, the memory requirements 
could become prohibitively large. Therefore, this approach was abandoned early in the present study. 

A convenient symbolic representation of the BB.BB scheme is obtained by differentiating Eq. (2.16) with 
respect to b to yield 

Q" n + 1 =Q" n -P"R" n -2P" , R' n -P""R“; n = 1, 2,3, ... (3.25) 

In contrast, the IIM for the SO Method (1) (DD.DD) can be represented compactly by combining Eqs. (3.11) and 
(3.12); the result is (by dropping the subscripts j,k) 

Q"m+1 _ Q«m _ pR"m, m= l i 2,3,... (3.26) 


Clearly, Eqs. (3.25) and (3.26) are symbolically equivalent only at convergence of the flow equations and 
the FO sensitivity equations. Computationally, however, Eqs. (3.25) and (3.26) are by no means equivalent; the 
potential for greater efficiency is with Eq. (3.26). This potential is seen from the previous discussion of the FO 
BB application of AD, where Eqs. (2.16) and (2.17) are compared. The potential for computational inefficiency is 
even greater now than previously for Eq. (2.16) because of the additional presence of the unwanted term P" n R n 
in Eq. (3.25) and because the unwanted term 2 P' n R' n represents a double computational evaluation of P' n R ,n . 

It is suggested that ADIFOR can be used to assist in the creation of the potentially efficient scheme represented 
by Eq. (3.26) in a manner that is completely analogous to the implementation of the AD-assisted FO method, 
ADII. Thus, the resulting method, known here as the ADII.SO scheme, involves the AD of all terms on the RHS 
only of Eq. (2.9), which is the residual R', of the linear FO sensitivity equations (Eq. (2.4)). These differentiated 
terms are then judiciously assembled for efficient operation on the RHS of the identical approximate LHS operator 
and algorithm that were also used to efficiently solve the nonlinear flow and the linear FO sensitivity equations. 

The ADH.SO scheme is potentially the most efficient implementation of the SO Method (1) (DD.DD) that 
is feasible, because HD construction of the scheme is too complex to be practical. It represents the only true 
SO IIM which, in principle, can be constructed in total with the present version of ADIFOR. (Recall that the SO 

* ADIFOR 2.0 has recently been installed and is being tested at NASA Langley with full Fortran support, improved error handling, and 
sparsity enhancements, but without SO SD capability. 
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Methods (2) and (3/4) require construction of the FO AV equations, in which the construction of some terms is 
not feasible via ADIFOR.) Actual implementation of the ADII.SO scheme is not included in the present study 
but is a topic of ongoing work. 

An alternative AD-assisted SO strategy is the BB application of ADIFOR to the scheme represented by Eq. 
(2.17); after differentiation with respect to b, the representation becomes 

Q"m+1 _ Q»m _ P B."m _ p' H.' m ; m=li2f 3,... (3.27) 

From the previous discussions of the analogous FO Eqs. (2.16) and (2.17), Eq. (3.27) clearly represents a significant 
improvement over Eq. (3.25) with respect to computational efficiency. However, Eq. (3.27) will also surely be 
less efficient than the true SO IIM (Eq. (3.26)). Symbolically (but not computationally), the two are equivalent 
only at convergence when R/ vanishes. 

The preceding method is constructed by the BB AD of source code which computes FO SD via an IIM (except 
the AVH scheme, for the present application). In principle, this code could be the source code for the previously 
discussed ADR scheme. The source code for the HD DDII method was selected in the sample problems, however, 
because it was more efficient. The SO method which results is known here as the DDII.BB scheme; by extension 
of this terminology, the FO DDIITC and DDIITCFR methods become the SO DDIITC.BB and DDIITCFR.BB 
schemes, respectively. By considering the latter two schemes, it is interesting to note that part of the original 
source code differentiated by ADIFOR is that which calculates the turbulence-modeling correction terms (i.e., 
those terms that were originally created by ADIFOR are then successfully differentiated by ADIFOR). This is the 
only example in the present work where this was actually attempted. 

For completeness, a final SO scheme is introduced here, called the AVH.BB method, which can be constructed 
via the BB application of ADIFOR to the source code for the FO AVH scheme. The AVn.BB scheme would yield 
inaccurate SO SD for the turbulent example problem because the full AD-generated correction of the turbulence- 
modeling terms could not be added to the FO AVII scheme. Symbolically, the AVH.BB scheme is represented 
by combining Eqs. (2.12) and (2.13) and differentiating with respect to b to yield 

A' m+1 = A' m - (P) T G' m - (P') T G m ; m= 1,2,3,... (3.28) 

In contrast, a similar representation of the HM for the SO Method (2) (AVX)D) is (Recall Eqs. (3.18) and (3.19)) 

A' m+1 = A' m - (P) T G' m ; m = 1, 2,3, ... (3-29) 
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Comparisons of the two schemes represented by Eqs. (3.28) and (3.29) yield conclusions with respect to 
computational efficiency that are analogous to those made previously for other FO and SO methods that employ 
BB applications of ADIFOR. 

The AVII.BB scheme, which essentially is an AD-assisted form of the SO Method (2) (AV DD), is not 
pursued further in this study because the SO Method (2) is always less efficient than the SO Method (3/4) 
(DDAV / A Vj\V) as demonstrated earlier. 

33 Computational Results: SO SD 

The FO SD results for the previous sample problems are extended here to include calculation of the complete 
Hessian matrices and These SO SD are calculated with different methods, and the results are 

compared on the basis of accuracy and computational time and memory. For the methods tested, considerable CW 
could have been saved by taking advantage computationally of the symmetry of the Hessian matrices; however, 
this was not done here in order to exploit this symmetry as an additional internal accuracy check. 

33.1 Accuracy Comparisons. The SO SD Hessian matrices are calculated for both the laminar and turbulent 
example problems with two basic methods known here as QA.CD and DDII.BB. The latter method has been 
described previously herein; more precisely specified, the DDII.BB scheme is applied to the laminar example, 
and the DDHTC.BB scheme is applied to the turbulent example. 

The QA.CD method can be described as a hybrid quasi-analytical/central finite-difference scheme, which is 
implemented in the following manner: 

1) Recall that for the FO CD method 12 machine-zero-converged solutions of the nonlinear flow equa- 
tions were generated: two perturbations for each design variable, a forward and a backward perturbation of 
Abj= ±5.0E-6 x bj. 

2) For each of the 12 perturbed nonlinear flow solutions, the complete set of FO SD were calculated using 
the FO HD QA code; specifically, the DDII and DDIITCFR methods were used for the laminar and turbulent 
examples, respectively. A reduction of 10 OM or better in the error of each linear system is specified; 72 linear 
systems are solved for these FO SD. 

3) The complete SO Hessian matrices were calculated with central finite-difference approximations using the 
QA FO derivatives calculated via that of the preceding discussion. 
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Thus, the QA.CD method for the SO SD performs the first differentiation exactly, and the second differentiation 
via a CD approximation. 

The SO SD that were calculated via the QA.CD method are presented in Tables (B.la) and (B.lb) (in 
Appendix B) for the laminar and turbulent sample problems, respectively. In these results, the actual numerical 
values of the derivatives are given for the main-diagonal and above-main-diagonal terms. The results presented 
for the below-main-diagonal terms are SD ratios, where each result shown has been divided by its equivalent 
above-main-diagonal term; clearly this calculation is an internal accuracy check of all the off-main-diagonal terms. 
As expected, these below-main-diagonal SD ratios are all unity to at least three and usually four or more significant 
digits. In these tables, the quasi-analytical first differentiation is with respect to bj (shown horizontally), followed 
by the CD approximate second differentiation with respect to bk (shown vertically). 

The SO results that were generated via the DDII.BB method are shown for the laminar example in Table 
(B.2a); similar results with the DDHTC.BB scheme for the turbulent example are shown in Table (B.2b). All of 
these results are given as SO SD ratios, where the numerical value of each result has been normalized by the 
numerical value of the respective term calculated via the QA.CD method. These tables clearly show the excellent 
agreement among the results obtained by these two methods, as expected. 

Prior to execution of the AD-enhanced code for the DDII.BB-type methods, the FO solution for the complete 
Q' matrix should be calculated first as input for the subsequent SO SD calculations. This initial Q' is calculated 
with the DDII-type methods (i.e., the original code from which the AD-enhanced code was created). An initial 
total of six linear systems is solved; a reduction of 10 OM in the error of each was specified. Subsequent execution 
of the AD-enhanced code produces the solutions of 36 linear systems for the complete SO SD Q"; an average 
reduction of 8 OM in the error of each system is specified in this case. These tight convergence tolerances should 
ensure accurate results, but at great expense in the computational timing comparisons presented subsequently. 

In addition to the 36 solutions for the SO terms, execution of the AD-enhanced code results in the solution 
of six linear systems for the FO Q'; this is a computationally wasteful, repeated solution for these terms, but 
their calculation cannot be avoided here because it is the function of the original code. The initial solution for 
Q' as input to the AD-enhanced code, discussed previously, can be avoided under certain conditions; however, 
such avoidance may not necessarily be efficient with respect to convergence rates and computer memory. This 
implementation is always possible and is also straightforward to invoke if the original code solves the FO sensitivity 
equations concurrently ; the initial input of the complete Q' matrix is replaced directly by the dynamic calculation 
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of these same terms, because they also evolve concurrently during the SO SD calculations. The original code used 
here solves the FO sensitivity equations sequentially, however. Consequently, avoidance of the initial Q' solution 
becomes a more complicated issue; in some cases where the complete Hessian matrices are not calculated, it is 
not even possible. An expanded discussion of the preceding ideas is planned for Ref. [20]. 

The complete set of terms on the main diagonal of the Hessian matrices can be calculated via the CD.CD 
method (i.e., the pure central finite-difference approximation of these terms) with the initial flow solution and 
the 12 perturbed flow solutions calculated previously for application in the FO CD and SO QA.CD schemes. 
An approximation of the complete set of cross-derivative terms in this manner is also possible, of course, but 
would require many additional perturbed flow solutions. Comparison of the results for these terms from the 
applications of the CD.CD and QA.CD methods is presented in Tables (B.3a) and (B.3b) for the laminar and 
turbulent examples, respectively; the results shown are SO SD ratios, where the CD.CD calculations reported here 
have been normalized by the respective QA.CD results. 

The two Tables, (B.3a) and (B.3b), show poor agreement among the results of these two methods. This dis- 
crepancy is attributed to inaccuracy in the CD.CD calculations; it occurs whenever a finite-difference perturbation 
is too small, the consequence of which is a necessity for accuracy in the function evaluations that is beyond the 
capability of the finite-arithmetic machine. With the present perturbation Abj= ±5.0E-6 x bj, this machine lim- 
itation was not exceeded for the accuracy desired in the FO CD results. For the SO CD.CD results, (Abj) 2 is now 
the divisor in the finite-difference expressions. Compared to the FO CD approximations, six or more additional 
significant digits of accuracy are required in the function evaluations to achieve the first digit of accuracy in the 
SO CD.CD results (with this perturbation). Of course, significant competing accuracy considerations emerge as 
the perturbation size is progressively increased. Larger perturbations were not tried for this demonstration because 
of the inordinate amount of CPU time required to produce SO derivatives, and the lack of any guarantee that the 
CD.CD results would be any better. Therefore, when the simple finite-difference method is applied, the selection 
of a good numerical perturbation size becomes progressively more difficult for higher order derivatives. 

33.2 Computational Time and Memory Comparisons. In this section, the SO methods of the previous 
section are further subdivided in a manner similar to that which was done previously for the various FO methods 
that were tested. As before, these subdivisions have no impact on the SD that are calculated but can greatly 
impact the computational time and memory requirements. 
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The QA.CD method is subdivided into two methods, depending on whether or not the FJ option is activated; 
this refers to the first phase of the method only, where the 12 nonlinear flow solutions are obtained. When this 
option is active, the method becomes the QA.CD(FJ) scheme; the method remains the simple QA.CD scheme 
when the option is inactive. 

The DDII.BB and DDIITC.BB methods, including the previously described DDHTCFR.BB efficiency op- 
tion that is associated with the latter method, are subdivided into the following six schemes: DDII.BB(6), 
DDn.BB(5+l), DDHTC.BB(6), DDIITC.BB(5+1), DDDTCFR.BB(6), and DDIITCFR.BB(5+1). These subdi- 
visions are made based on different applications of the parameter g$p$ in the AD-enhanced code, as described 
previously for the FO results. 

Comparisons of total CPU times are shown in Table (B.4) for the laminar and turbulent sample problems. 
All reported timings do not include the cost of the initial flow solution. For the QA.CD and QA.CD(FJ) schemes, 
the cost for the 12 perturbed nonlinear flow solutions and the 72 linear system derivative solutions is included in 
the reported timings. For the AD-assisted DDII.BB-type schemes, the computational cost is included for initially 
solving the FO equations for all Q'. (Recall that this FO solution is required input at the outset of the execution 
of the AD-enhanced code for the SO SD.) Table (B.5) shows comparisons (excluding the QA.CD and QA.CD(FJ) 
methods) of CPU time per iteration per SO linear system solved. The cost of the initial FO solution for Q' is 
not included in the results of this table. As described previously for the FO results, the individual effects of g$p$ 
= 5 and g$p$ = 1 are shown in Table (B.5). Finally, the total computer memory requirements for the different 
methods are compared in Table (B.6). 

3.4 Conclusions: SO SD 

Conclusions based on the calculations for the SO SD are enumerated subsequently. 

1) The calculation of SO aerodynamic SD by a pure finite-difference scheme is much more sensitive to the 
selection of a proper perturbation size than is the calculation of the FO SD; this difficulty is noted in addition to the 
extreme computational cost of the method, particularly if the complete Hessian matrices (with cross-derivatives) 
are computed. These difficulties for the SO derivatives can be mitigated to a large extent if the FO derivatives, 
at least, are calculated via one of the HD or AD methods discussed previously. 

2) The good agreement that is seen in the AD-assisted SO results (when compared with the results from the 
QA.CD method) confirms that ADIFOR can be successfully used as a tool to construct and implement the SO 
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methods from source code that evaluates the FO derivatives. In principle, the source code for FO derivatives 
can be constructed either in whole or in part by either HD or AD. Without ADIFOR, the SO schemes could not 
be implemented even for simple CFD codes; construction of the source code by hand would not be feasible to 
evaluate the extremely large number of complex SO terms that are involved. 

3) The previously discussed limitations of ADIFOR when the FO AVH methods are constructed carries over at 
least indirectly in all the SO methods except Method (1) (DDDD) because the other three SO methods require 
the FO AV equations to be solved. This requirement is significant because Method (3/4) is potentially by far the 
most efficient method for particular combinations of NDV and NOF. 

4) As expected, the computational results, although accurate, were very costly in terms of CPU time and (for 
the DDII.BB-type methods) computer memory. The possibility exists that these computational costs can be 
significantly reduced, however, as described in the subsequent conclusions; this possibility is also a topic of 
ongoing study. 

5) For the particular combination of NDV and NOF studied here. Method (3/4) is potentially far more efficient 
than the methods actually tested; it is projected to cost only slightly more than the cost of solving the large 
systems of equations for the FO DD and AV schemes — a total of only nine large systems. This cost, taken from 
the HD FO results, was a total of only 583 and 1014 CPU seconds for the laminar and turbulent-flow examples, 
respectively. The results for the turbulent-flow example would have inaccuracies traced to the FO AVII scheme 
(which could not be corrected via AD for all of the turbulence-modeling terms). 

6) In contrast with the nine large system solutions that would be required for Method (3/4), the DDflLBB type of 
implementation of Method (1) included the solution of 42 large systems — six for FO derivatives and 36 relatively 
inefficient (as discussed subsequently) solutions for the SO derivatives. Note that by simply taking advantage 
computationally of the symmetry of the Hessian matrices, the number of solutions for SO derivatives could be 
reduced from 36 to 21. 

7) As discussed previously, the DDII.BB-type approach is not the most efficient implementation of Method (1) 
because it is not a true IIM; it is believed that the true SO IIM implementation, known here as the ADII.SO 
approach, would be the most efficient feasible implementation of Method (1). The improvement in CPU time and 
memory requirements that results from the ADII.SO implementation (compared to the DDII.BB-type approach) 
is projected to be roughly proportional to those improvements seen in a comparison of the two FO methods ADII 
and BBFJ (that is, ADII.SO could be possibly 40%-70% faster than DDII.BB). 
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8) The computational efficiency, in terms of CPU time and memory, is affected greatly by the selection of the 
parameter g$p$ in the AD-assisted SO schemes. This effect is similar to that noted previously in the FO results. 


4.0 Summary and Final Conclusions 

A number of different approaches for computing first-order (FO) aerodynamic sensitivity derivatives (SD) 
have been tested, and the results have been compared on the basis of accuracy, computational time, and computer 
memory requirements. The methods represent a broad spectrum of choices: finite-difference methods, hand- 
differentiated (HD) incremental iterative method (EM) schemes, and black-box (BB) automatic differentiation 
(AD) methods. In addition, combinations and variations of these are possible. 

The automatically differentiated incremental iterative (ADII) scheme stands out as a very well-rounded 
combination of the best features of the other methods. Although it is not yet as efficient as the HD schemes, 
it is more efficient than all other methods tested — substantially more efficient than the simple BB and central 
finite-difference (CD) approaches. The ADII method is not quite as easy to implement as the BB approach. 
However, when compared with the HD methods, the ADII scheme can be implemented quickly, reliably, and 
with very accurate results, even for very complex computational fluid dynamics (CFD) codes. The ADII scheme 
also requires less expansion of computer memory than some of the other AD-assisted methods. When compared 
with the best CD method, the most efficient implementation of the ADII scheme improved computational efficiency 
by a factor of approximately 3.6 for the laminar example, and 1.7 for the turbulent example; the memory increase 
was about 1.6 times that of the original flow code. The computational penalty which is associated with the 
AD-assisted differentiation of the turbulence-modeling terms is significant and disproportionately large for all 
applicable methods; this requires additional examination for possible improvements. 

A complete procedure has been described by which the discrete second-order (SO) aerodynamic SD can 
be calculated from modern CFD codes. The assistance of AD is required for the implementation, because of 
the extremely large number of complex terms that are required. Initially, four SO methods were presented; in 
effect, this selection was reduced to only two, because SO Methods (3) and (4) were shown to be equivalent and 
unconditionally more efficient than Method (2). Thereafter, the choice between either Method (1) or Method (3/4) 
is dependent upon the relative size of NDV and NOF; a large NOF favors Method (1). The specific criterion to 
be used for this selection was defined herein. 
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The SO Method (3/4) requires that the large systems for the FO DD and AV schemes are solved first. 
Thereafter, the computation of the SO SD is noniterative (i.e., additional simultaneous solutions of large systems 
are not needed). Method (1) requires that the large systems for the FO DD scheme are solved first; the FO AV 
systems are not solved, which can be a significant advantage over Method (3/4) for reasons that have been detailed 
herein. Thereafter, the SO SD are computed by solving large systems for SO terms; however, these equations 
and the efficient solution methods for them are analogous to those for the FO DD method. The computational 
results for the SO SD were highly accurate, which confirms that ADIFOR is a reliable tool in constructing these 
SO methods. Significantly improved results with respect to computational efficiency for SO SD are expected in 
the future from ongoing work. 
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7.0 Appendix A: Tables for First-Order Results 


Table (A.l): First-Order Sensitivity Derivatives and Sensitivity-Derivative Ratios; Laminar Example 



Solution 

Design Variable 

dC L 

dC D 

dCu 


Method 

bj 

db 

db 

db 

— 

Central Differences 

T 

-1.392E+00 

+2.019E-01 

+1.805E-01 


Method (CD) 

c 

-H5.583E+00 

+7.583E-02 

-2.240E+00 

— 

SDcd 

L 

-1.154E-02 

+5.540E-05 

-2.122E-02 



a 

+6.122E+00 

+9.181E-02 

-3.166E-02 



Moo 

+5.438E-03 

+1.628E-02 

-4.732E-03 



Re 

+5.958E-06 

-4.912E-06 

-6.563E-07 


Direct 

T 

1.0000 

1.0000 

1.0000 


Differentiation 

Incremental 

c 

1.0000 

1.0000 

1.0000 


Iterative Method 

L 

1.0000 

1.0009 

1.0000 

— 

(DDII) 

a 

1.0000 

1.0000 

1.0004 


SP DDlt 

SDcd 

Moo 

1.0000 

1.0000 

1.0000 



Re 

1.0000 

1.0000 

1.0000 


Adjoint- Variable 

T 

1.0000 

1.0000 

1.0000 


Incremental 

c 

1.0000 

1.0000 

1.0000 


Iterative Method 





(AV1I) 

L 

1.0000 

1.0009 

1.0000 


SD*vn 

a 

1.0000 

1.0000 

1.0004 


SDcd 

Moo 

1.0000 

1.0000 

1.0000 



Re 

1.0000 

1.0000 

1.0000 

— 

Automatic 

T 

1.0000 

1.0000 

1.0000 


Differentiation in 
Incremental 

C 

1.0000 

1.0000 

1.0000 

— 

Iterative Form 

L 

1.0000 

1.0009 

1.0000 


(ADII) 

a 

1.0000 

1.0000 

1.0004 

— 

SP AD1I 

SDcd 

Moo 

1.0000 

1.0000 

1.0000 



Re 

1.0000 

1.0000 

1.0000 


Automatic 

T 

1.0000 

1.0000 

1.0000 


Differentiation, 

c 

1.0000 

1.0000 

1.0000 


’’Black Box” 





Method (BB) 

L 

1.0000 

1.0009 

1.0000 


SDrb 

a 

1.0000 

1.0000 

1.0004 


SDcd 

Moo 

1.0001 

1.0000 

1.0000 

— 


Re 

1.0000 

1.0000 

1.0000 
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Table (A .2): First-Order Sensitivity Derivatives and Sensitivity-Derivative Ratios; Turbulent Example 


Solution 

Design Variable 

dC L 

dC D 

dC M 

Method 

bj 

db 

db 

db 

Central Differences 

T 

+7.919E-01 

+2.744E-01 

-4.153E-01 

Method (CD) 

c 

+2.063E+01 

+6.776E-01 

-5.770E+00 

SDcd 

L 

+1.108E-01 

-1.174E-02 

-5.350E-02 


a 

+1.300E+01 

+4.346E-01 

-6.328E-01 


m to 

+2.040E+00 

+1.969E-01 

-5.972E-01 


Re 

-1.185E-09 

-2.829E-10 

+1.497E-10 

Direct 

T 

0.2874 

0.9672 

0.7523 

Differentiation and 
Adjoint-Variable 

c 

0.9415 

0.9609 

0.9560 

Incremental 

L 

12077 

0.9806 

1.0447 

Iterative Methods, 
Uncorrected (DDII 

a 

0.9221 

0.9663 

0.7388 

and AVII) 

SDnpn , nr < SP AY1 j 
SDcd ^ SD C d 

Moo 

0.8690 

0.9754 

1.7251 

0.9093 

-2.9367 

Re 

-3.4966 

Direct 

T 

1.0000 

1.0000 

1.0000 

Differentiation 

Incremental 

c 

1.0000 

1.0000 

1.0000 

Iterative Method, 

L 

1.0000 

1.0000 

1.0000 

Corrected 

(DDIITC) 

a 

1.0000 

1.0000 

1.0000 


Moo 

1.0000 

1.0000 

1.0000 

SDddiitc 




1.0000 

SDcd 

Re 

1.0000 

1.0000 

Automatic 

T 

1.0000 

1.0000 

1.0000 

Differentiation in 
Incremental 

C 

1.0000 

1.0000 

1.0000 

Iterative Form 

L 

1.0000 

1.0000 

1.0000 

(ADIT) 

a 

1.0000 

1.0000 

1.0000 

SDadii 

SD C d 

Moo 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 

Automatic 

T 

1.0000 

1.0000 

1.0000 

Differentiation, 
"Black Box" 

C 

1.0000 

1.0000 

1.0000 

Method (BB) 

L 

1.0000 

1.0000 

1.0000 

SDbb 

a 

1.0000 

1.0000 

1.0000 

SDcd 

Moo 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 
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Table (A .3): Computational Timing Comparisons: Total CPU Time (secs) 


Method 

Laminar 

Turbulent 

CD 

9,800* 

12,000* 

CDFJ 

2,962 

4,057 

AVII 

324 

300 

DDE 

259 

590 

DDHTC 

Not Applicable 

2,000* 

DDIITCFR 

Not Applicable 

714 

ADII(6) 

1,722 

4,879 

ADn(5+l) 

820* 

2,300* 

BB(6) 

14,000* 

27,000* 

BB(5+1) 

12,000* 

24,000* 

BBFJ(6) 

2,696 

6,964 

BBFJ(5+1) 

1,400* 

3,600* 


"Projected result based on timings from Table (A.4) 


Table (A.4): Computational Timing Comparisons: CPU Time (secs) / Iteration / Linear System Solved 


Method 

Laminar 

Turbulent 

AVII 

0.06196 

0.06293 

DDII 

0.06325 

0.06314 

DDHTC 

Not Applicable 

0.2285 

DDIITCFR 

Not Applicable 

0.07996 

ADII(6) 

0.4196 

0.5432 

ADII(5) \ ADII(l) 

0.1705 \ 0.3453 

0.2096 \ 0.5107 

BB(6) 

3.448 

3.570 

BB(5)\BB(1) 

2.033 \ 8.010 

2.050 \ 8.096 

BBFJ(6) 

0.7694 

0.9046 

BBFJ(5) \ BBFJ(l) 

0.3039 \ 0.8827 

0.3469 \ 1.046 
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Table (A.5): Total Computer Memory Comparisons 


Method 

Total Memory (Mw) 

CD, CDFJ 

5.27 

DDII, AVTI, DDHTC, DDHTCFR 

7.39 

ADII(6) 

9.07 

ADII(5) \ ADH(l) 

8.36 \ 5.09 

BB(6), BBFJ(6) 

34.65 

BB(5), BBFJ(5) \ BB(1), BBFJ(l) 

29.94 \ 10.27 
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8.0 Appendix B: Tables for Second-Order Results 


Table (B.la): Second-Order Sensitivity Derivatives and Sensitivity-Derivative Ratios: QA.CD Method; 


Laminar Example 



b k \bj 

T 

c 

L 

a 

Moo 

Re 


T 

-2.280E+01 

-8.390E+01 

-6.579E-01 

-3.661E+01 

-4.728E+00 

+3.752E-06 


c 

1.0000 

-6.708E+01 

+1.029E-01 

-3.048E401 

-5.473E+00 

+3.680E-04 

dbkdbj 


L 

1.0000 

1.0000 

-1.079E-02 

+1.803E-01 

-1.351E-01 

+4.409E-06 

\SD 

a 

1.0000 

1.0000 

1.0000 

-2.376E+01 

+1.700E+00 

+1.043E-04 

SD\ 

Ratios 

Moo 

1.0000 

1.0000 

1.0000 

1.0000 

-1.222E-01 

+7.902E-06 


Re 

1.0001 

1.0000 

1.0000 

1.0000 

1.0000 

-1.446E-09 


T 

+1.180E+00 

+ 1.017E+00 

-2.651E-02 

-4.182E-01 

+1.690E-01 

-5.236E-06 


C 

1.0000 

+9.983E+00 

-6.234E-02 

-1.352E+00 

+2.412E-01 

-2.481E-06 

d Cd 
dbkdbj 


L 

1.0000 

1.0000 

+4.952E-03 

+3.478E-02 

-3.804E-03 

+6.470E-08 

\SD 

a 

1.0000 

1.0000 

1.0001 

+5.606E+00 

+9.573E-02 

+2.326E-06 

SD\ 

Ratios 

Moo 

1.0000 

1.0000 

1.0000 

1.0001 

+4.447E-02 

-3.945E-07 


Re 

1.0000 

1.0001 

1.0000 

1.0000 

1.0000 

+1.606E-09 


T 

+2.965E+00 

+1.269E+01 

+1.907E-01 

+3.047E+00 

+6.144E-01 

+4.748E-07 


c 

1.0000 

+8.244E+00 

+1.233E-02 

+3.545E+00 

-6.423E-01 

-5.805E-05 

db k dbj 


L 

1.0000 

1.0000 

-1.076E-03 

+3.947E-02 

+1.459E-02 

-1.481E-06 

\SD 

a 

1.0000 

1.0000 

1.0000 

+2.206E+00 

+2.385E-01 

+6.425E-07 

SD\ 

Ratios 

Moo 

1.0000 

1.0000 

1.0000 

0.9999 

-3.726E-02 

-2.892E-07 


Re 

1.0000 

1.0000 

1.0000 

0.9999 

0.9998 

+1.547E-10 
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Table (B.lb): Second-Order Sensitivity Derivatives and Sensitivity-Derivative Ratios: QA.CD Method; 


Turbulent Example 



bk\bj 

T 

c 

L 

a 

Moc 

Re 


T 

-1.587E+02 

-5.480E+02 

-5.168E-01 

-1.883E+02 

-1.266E+02 

-1.664E-09 


C 

1.0000 

-1.780E+03 

-2.415E+00 

-6.383E+02 

-4.057E+02 

-3.025E-08 

db k dfcj 


L 

1.0000 

1.0000 

-5.723E-02 

+3.011E-01 

-7.414E-01 

-1.044E-09 

\SD 

a 

1.0000 

1.0000 

1.0000 

-1.772E+02 

-1.068E+02 

-1.355E-08 

SD\ 

Ratios 

Moo 

1.0000 

1.0000 

1.0000 

1.0000 

-1.020E+02 

+6.374E-09 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

+1.512E-16 


T 

+7.245E+00 

+1.763E+01 

-2.782E-01 

+1.278E+01 

+6.118E+00 

-7.085E-10 


c 

1.0000 

+6.397E+01 

-6.362E-01 

+3.285E+01 

+1.730E+01 

-3.797E-09 

db k dbj 

L 

1.0000 

1.0000 

+1.683E-02 

-5.967E-01 

-2.184E-01 

+1.675E-11 

\SD 

a 

1.0000 

1.0000 

1.0000 

+2.879E+01 

+1.235E+01 

-1.916E-09 

SD\ 

Ratios 

Moo 

1.0000 

1.0000 

1.0000 

1.0000 

+5.029E+00 

-1.054E-09 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

+1.104E-16 


T 

+4.993E+00 

+ 1.492E+01 

+3.707E-01 

-1.283E+01 

-5.695E+00 

+3.513E-09 


C 

1.0000 

-3.965E+00 

+8.127E-01 

-4.331E+01 

-2.995E+01 

+1.753E-08 



L 

1.0000 

1.0000 

-3.088E-03 

+1.308E+00 

+3.977E-01 

+2.034E-10 

VSD 

a 

1.0000 

1.0000 

1.0000 

-3.227E+01 

-2.348E+01 

+4.208E-09 

SD\ 

Ratios 

Moo 

1.0000 

1.0000 

1.0000 

1.0000 

-1.173E+01 

+2.259E-09 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

-1.668E-17 
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Table (B.2a): Second-Order Sensitivity-Derivative Ratios: DDILBB Method, Normalized by Results 


from the QA.CD Method; Laminar Example 



b k \bj 

T 

c 

L 

a 

Moo 

Re 


T 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


C 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

d a C L 

db k dbi 


L 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

so 

SD 

a 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

Ratios 


Moo 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


T 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


c 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

d Cd 
dbkdDt 


L 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

so 

SD 

a 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0000 

Ratios 


Moo 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


T 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


c 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

&7 


L 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

SO 

SD 

a 

1.0000 

1.0000 

1.0000 

1.0000 

0.9999 

0.9999 

Ratios 


M « 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0002 

1.0000 
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Table (B.2b): Second-Order Sensitivity-Derivative Ratios: DDIUC.BB Method, Normalized by Results 
from the QA.CD Method; Turbulent Example 



b k \bj 

T 

c 

L 

a 

Moo 

Re 


T 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


c 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

d*Ci. 

dbirdbj 

L 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 

so 

SD 

a 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

Ratios 


Moo 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


T 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


c 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

d a Cp 

dbk dbi 


L 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

so 

SD 

a 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

Ratios 


Moo 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


T 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


c 

1.0000 

0.9999 

1.0000 

1.0000 

1.0000 

1.0000 

L 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 

so 

SD 

a 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

Ratios 


Moo 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


Re 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 
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Table (B3a): Second-Order Sensitivity-Derivative Ratios; CD.CD Method Normalized by Results from 
the QA.CD Method; Main-Diagonal Terms Only; Laminar Example 


— 


T 

C 

L 

a 

Moo 

Re 


d 3 C L 
db f 

-0.7096 

-18.31 

- 71.90 

-18.78 

+2.138 

-3.023 

— 

d 3 Cp 
db Y 

+1.803 

44.341 

+5.158 

+3.027 

+1.214 

+1.095 


<J j Cm 

dbj 

+1.017 

+5.160 

-18.00 

+0.6872 

+0.9608 

+0.6613 

— 

Table (B3b): Second-Order Sensitivity-Derivative Ratios; CD.CD Method Normalized by Results from 
the QA.CD Method; Main-Diagonal Terms Only; Turbulent Example 



T 

C 

L 

a 

Moo 

Re 


d 3 C L 

dbj 

+0.8072 

+0.5389 

-12.55 

-1.274 

+1.000 

+31.81 

— 

d 3 Cp 

dbj* 

+1.028 

+0.8928 

+1.390 

+1.072 

+0.9992 

+1.353 

— 

d a CM 

dbj 3 

+1.569 

-35.16 

-12.87 

+0.2042 

+0.9966 

-17.68 
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Table (B.4): Computational Timing Comparisons: Total CPU Time (secs) 


Method 

Laminar 

Turbulent 

QA.CD 

14,000* 

23,000* 

QA.CD(FJ) 

7,018 

14,941 

DDII.BB(6) 

14,543 

Not Tested 

DDH.BB(5+1) 

4,300* 

Not Tested 

DDIITC.BB(6) 

Not Applicable 

64,000* 

DDIITC.BB(5+1) 

Not Applicable 

51,000* 

DDHTCFR.BB(6) 

Not Applicable 

33,899 

DDIITCFR.BB(5+1) 

Not Applicable 

14,000* 


‘Projected result based on timings from Table (B.5) 


Table (B.5): Computational Timing Comparisons: CPU Time (secs) / Iteration / Linear System Solved 


Method 

Laminar 

Turbulent 

DDII.BB(6) 

0.57562 

Not Tested 

DDII.BB(5) \ (1) 

0.1409 \ 0.2675 

Not Tested 

DDIITC.BB(6) 

Not Applicable 

1.2295 

DDIITC.BB(5) \ (1) 

Not Applicable 

0.6630 \ 2.5500 

DDIITCFR.BB(6) 

Not Applicable 

0.6474 

DDIITCFR.BB(5) \ (1) 

Not Applicable 

0.2016 \ 0.5216 


Table (B.6): Total Computer Memory Comparisons 


Method 

Total Memory (Mw) 

QA.CD, QA.CD(FJ) 

7.39 

DDII.BB(6), DDIITC.BB(6), DDIITCFR.BB(6) 

48.85 

DDn.BB(5) \ (1), DDHTC.BB(5) \ (1), 

42.38 \ 14.55 

DDIITCFR.BB(5) \ (1) 
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9.0 Appendix C: Tables of Acronyms 


Table (C.1): General Terms 



2-D 

Two-dimensional 


3-D 

Three-dimensional 


AD 

Automatic differentiation 

— 

ADIFOR 

Automatic DIFferentiation of FORtran (software tool) 


CFD 

Computational fluid dynamics 

— 

CPU 

Central processing unit (of computer) 


CW 

Computational work (CPU time and memory used) 

_ 

FJ 

Frozen Jacobian (option for fast convergence of a CFD solution) 


FO 

First-order (derivatives or Taylor series expansions) 


FR 

Indicates that the turbulence correction terms are frozen 


g$p$ 

Parameter in all ADIFOR -generated FORTRAN source code 
governing the length of do-loops to calculate a set of derivatives 

— 

HD 

Hand-differentiated (as opposed to finite differences or automatic 
differentiation) 

— 

IIM 

Incremental Iterative Method (as opposed to standard methods, 
for solution of linear matrix equations) 


LHS 

Left-hand side (of an equation) 

— 

MDO 

Multidisciplinary design optimization 


NDV 

Number of design variables 



NOF 

Number of output functions 


OM 

Order of magnitude 

— 

QA 

Quasi-analytical (differentiation of discretized equations by 
calculus, or '’analytically", and by hand) 


RHS 

Right-hand side (of an equation) 

— 

SD 

Sensitivity derivative(s) 


SO 

Second-order (derivatives or Taylor series expansions) 

— 

TLNS 

Thin-layer Navier-Stokes (equations or code) 


TC 

Indicates that AD-generated turbulence-correction terms are 
included 
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Table (C*2): First-Order Formulations and Solution Strategies 


ADII 

Automatic differentiation in incremental iterative form 

AV 

Adjoint-variable (formulation for derivative equations) 

AVII 

Adjoint-variable incremental iterative (strategy) 

aviitc 

AVH with AD-generated turbulence-modeling correction (not 
currently possible) 

BB 

Black-Box (straightforward application of AD, as opposed to an 
application in incremental iterative form) 

BBFJ 

Black-Box with a frozen Jacobian 

CD 

Central finite differences 

CDFJ 

Central finite differences with a frozen Jacobian 

DD 

Direct differentiation (as opposed to adjoint-variable formulation) 

DDII 

Direct differentiation incremental iterative (strategy) 

DDIITCFR 

DDIITC with the AD-generated turbulence-modeling correction 
frozen 

DDIITC 

DDII with AD-generated turbulence-modeling correction 

*(6) 

Implementation for six design variables concurrently 

*(5+1) 

Implementation first for five design variables concurrently, 
followed by one design variable 


* Indicates each of the methods ADD, BB, and BBFJ 
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Table (C3): Second-Order Formulations and Solution Strategies 



^ADII.SO 

SO automatic differentiation in incremental iterative form 


AVDD 

FO adjoint-variable approach, then SO direct differentiation 


AVAV 

FO adjoint-variable approach, then SO adjoint-variable approach 

— 

Uvn-BB 

FO AVn, then black-box automatic differentiation 


tBB.BB 

Black-box automatic differentiation applied twice 

— 

CD.CD 

Second-order central finite differences 


DDAV 

FO direct differentiation, then SO adjoint-variable approach 


DDDD 

FO direct differentiation, then SO direct differentiation 


DDII.BB 

FO DDn, then black-box automatic differentiation 


DDIITC.BB 

FO DDEITC, then black-box automatic differentiation 


DDIITCFR.BB 

FO DDHTCFR, then black-box automatic differentiation 


QA.CD 

FO quasi-analytical differentiation, then first-order central finite 
differences 


QA.CD(FJ) 

QA.CD with a frozen Jacobian 


*(6) 

Implementation for six design variables concurrently 


*(5+1) 

Implementation first for five design variables concurrently, 
followed by one design variable 


* Not implemented. 

* Indicates each of the methods DDII.BB, DDIITC.BB, and DDIITCFR.BB 


50 


